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Abstract 

A  non-isothermal  and  three-dimensional  numerical  model  of  a  PEM  fuel  cell  was  developed  to  compute  the  water  and  heat  management. 
Transport  of  water  in  the  polymer  membrane,  phase  change  of  water  in  the  cathode  porous  medium  and  capillary  flow  to  the  gas  channels  were 
determined.  Influences  of  these  phenomena  on  fuel  cells  and  conditions  that  may  affect  their  performance  have  been  numerically  evaluated.  Output 
variables  are  velocity,  temperature,  mass  fraction,  current  density,  voltage  loss,  water  content  of  the  polymer  membrane,  saturation  and  liquid 
flow  fields.  Cell  voltage  and  total  current  density  of  PEM  fuel  cell  were  computed  as  well.  Results  show  that  there  may  be  severe  mass  transfer 
limitations  depending  either  on  the  design  or  on  the  water  management  of  the  cell.  For  the  chosen  conditions,  the  polymer  membrane  can  keep 
and  even  increase  its  water  content,  as  long  as  inlet  flows  are  injected  at  100%  relative  humidity.  In  case  the  fuel  cell  is  operated  under  dehydrating 
conditions,  the  decrease  of  the  water  content  of  the  polymer  electrolyte  may  affect  the  performance.  The  variations  of  temperature  were  small. 
However,  temperature  plays  an  important  role  in  the  cathode  reaction  rate  of  the  cell  and  in  the  dehydration  of  the  polymer  membrane.  Numerical 
results  and  experimental  data  were  found  to  be  in  good  agreement. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

The  performance  of  PEM  fuel  cells  depends  highly  on  the 
water  content  of  the  polymer  membrane,  which  is  inversely  pro¬ 
portional  to  the  proton  resistance.  Therefore,  decreasing  water 
content  leads  to  a  diminishing  of  performance  due  to  the  increase 
of  the  voltage  loss  in  the  membrane.  In  addition  to  this  point, 
dehydrating  conditions  in  PEM  fuel  cells  may  not  only  affect  the 
proton  resistance  of  the  membrane,  but  also  the  proton  resistance 
of  the  polymer  electrolyte  contained  in  the  cathode  catalyst  layer. 
Then,  the  cathode  reaction  rate  is  affected  by  changes  of  local 
cathode  overpotential  caused  by  high  proton  resistance  [1,2]. 

Therefore,  to  avoid  dehydration  of  the  polymer  electrolyte, 
the  inlet  flows  are  usually  injected  as  a  mixture  with  water  vapor 
at  a  certain  temperature.  However,  produced  water  and  lack  of 
diffusion  flux  for  the  disposal  of  this  water  may  cause  conden¬ 
sation  in  the  cathode  porous  medium,  which  may  clog  pores  and 
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block  the  transport  of  species  to  reactive  sites.  For  these  reasons, 
the  management  of  water  in  the  PEM  fuel  cell  becomes  highly 
important. 

In  order  to  achieve  a  complete  solution  of  the  transport  phe¬ 
nomena  in  PEM  fuel  cells,  the  water  management,  which  in  turn 
requires  the  consideration  of  the  heat  management,  must  be  ac¬ 
curately  taken  into  account.  In  order  to  establish  the  influence  of 
the  water  on  the  performance  of  the  fuel  cell,  the  water  move¬ 
ment  in  the  polymer  membrane  as  well  as  the  condensation  and 
capillary  flow  in  the  cathode  gas  diffusion  layer  and  catalyst 
layer  should  be  properly  estimated.  Finally,  the  flow  field  in  the 
gas  channel  should  be  able  to  dispose  the  condensed  water.  If 
these  aspects  are  neglected,  numerical  results  may  be  erroneous, 
given  that  the  behavior  of  the  whole  system  is  highly  dependant 
on  the  liquid  and  gaseous  water  transport. 

Many  approaches  for  calculation  of  one,  two  and  three  di¬ 
mensional  variations  of  important  properties  inside  PEM  fuel 
cells  have  been  made.  However,  each  approach  put  emphasis  on 
different  aspects.  Wang’s  review  [3]  shows  a  complete  descrip¬ 
tion  of  the  most  important  models  created  for  the  calculation  of 
different  variables  involved  in  PEM  fuel  cells.  Up  to  now,  some 
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Nomenclature 

Alg/Vl 

gas-liquid  interfacial  area  (m-1) 

^superf 

Pt  area  per  Pt  mass  (m2  kg-1) 

ARH 

anode  relative  humidity 

cp 

volumetric  thermal  capacity  (Jkg-1  K-1) 

C 

molar  concentration  (kmolm-3) 

^C^ref 

oxygen  reference  molar  concentration 

(kmol  m-3) 

Q)2S 

oxygen  molar  concentration  at  the  agglomerate 
surface  (kmolm-3) 

CRH 

cathode  relative  humidity 

D 

diffusivity  (m2  s-1) 

F 

Faraday  constant 

8 

gravity  acceleration  (ms-2) 

Gref 

reference  exchange  current  density  (Am-2) 

k 

thermal  conductivity  (W  m-1  K-1) 

Kp 

hydraulic  permeability  (m2) 

m^i 

Pt  mass  per  volume  (kg  m-3) 

M 

molecular  weight  (kg  kmol-1) 

^glh2o 

volumetric  condensation  rate  (kg  s-1m-3) 

N 

molar  flux  (mol  m-2  s-1) 

^Drag 

electro-osmotic  drag  factor  (molH2o(molH+)_1) 

P 

pressure  (Pa) 

r 

reaction  rate  (kmol  s-1  m-3) 

R 

universal  constant  of  gases 

^agg 

agglomerate  ratio  (m) 

liquid  saturation 

Si 

source  term  of  species  i 

Sm 

source  term  of  mass 

Sp 

source  term  of  energy 

Sy 

source  term  of  momentum 

T 

temperature  (K) 

V 

velocity  (ms-1) 

X 

molar  fraction 

y 

mass  fraction 

Greek 

ot  c 

transfer  coefficient  cathode  side 

€ 

porosity 

0 

cathode  overpotential  (V) 

0 

contact  angle  (rad) 

X 

water  content  (molH2o(fnolso-)-1) 

M 

dynamic  viscosity  (kgm-1  s-1) 

Q 

density  (kgm-3) 

a 

surface  tension  (N  m-1) 

Subscripts 

agg 

agglomerate 

c 

capillary 

G 

gas 

H+ 

proton 

H20 

water 

i 

species  i 

J 

species  j 

L 

liquid 

m 

mixture 

o2 

oxygen 

s 

solid 

so3- 

sulfonate  group 

models  have  been  provided  to  estimate  the  movement  of  liquid 
and  gaseous  water  in  the  cathode  porous  medium  of  the  fuel  cell 
[2,4-7]. 

The  objective  of  this  work  is  to  present  a  numerical  solution  of 
the  water  and  heat  management  in  PEM  fuel  cells.  The  polymer 
membrane  is  considered  as  a  component  highly  affected  by  its 
external  conditions,  hence  the  mass  transfer  between  the  mem¬ 
brane,  the  cathode  and  the  anode  porous  mediums  is  taken  into 
account,  in  order  to  determine  the  actual  water  content  field,  so 
the  actual  voltage  losses  and  the  cathode  electrolyte  overpoten¬ 
tial  are  obtained.  The  liquid  water  saturation  factor  in  the  porous 
medium  is  estimated  by  relations  based  on  the  kinetic  theory  of 
gases  and  by  the  Leverett  function  for  the  capillary  movement 
of  liquid  in  unsaturated  hydrophobic  porous  mediums  [2, 4, 6, 7]. 

2.  Model  description 

A  simple  straight  PEM  fuel  cell  was  taken  as  computational 
domain,  given  that  the  main  purpose  of  this  work  is  to  study 
the  transport  behavior  of  the  species  in  each  component.  Using 
a  symmetry  assumption,  the  domain  was  vertically  split  in  two 
parts  and  only  the  left  part  was  calculated  (Fig.  1).  The  model 
is  three-dimensional  and  simultaneously  calculates  the  behavior 
of  gas  channels,  catalyst  layers  and  polymer  membrane. 

2.7.  Fundamental  equations 

In  order  to  obtain  the  velocity  fields,  the  Navier-Stokes  equa¬ 
tions  were  used  in  the  following  form: 

V  •  VqV  =  -VP  +  V/iVV  +  SV  (1) 


along-channel 


Fig.  1.  Control  volume  for  the  numerical  solution. 
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The  energy  balance  is  expressed  as 

V  •  VqcpT  =  -VkVT  +  ST  (2) 

and  the  species  mass  balance  is  represented  by 

V  •  Vgyi  =  -VDUV^  +  Si  (3) 

Finally,  each  equation  is  coupled  with  the  continuity  equa¬ 
tion: 

=  Sm  (4) 


inverse  of  the  condensation  process.  The  gas-liquid  interfacial 
area  Alg/  Vl  was  arbitrarily  assumed  as  constant. 

All  equations  of  balance  are  discretisized  and  solved  simul¬ 
taneously  using  the  Jacobi-Gauss-Seidel  relaxation  iterative 
numerical  method.  This  method  and  boundary  conditions  were 
written  in  Fortran90  for  calculation  of  the  output  parameters  of 
the  model. 

2.2.  Transport  properties 


For  the  polymer  membrane,  a  water  content  balance  was 
made  in  order  to  consider  electro-osmotic  drag  and  back  dif¬ 
fusion  of  water  molecules: 

MSO-  dA'oraaVr 

VDxVX  = - 5 - —  (5) 

Bso-  dx 


Two-phase  flow  in  the  cathode  porous  medium  (gas  diffusion 
layer  and  catalyst  layer)  is  assumed  as  separated  phases.  Berning 
and  Djilali  [6]  proposed  a  liquid  water  velocity  field  considering 
both  convection  and  capillary  forces.  In  this  work,  the  liquid  is 
assumed  to  flow  individually  driven  by  capillary  pressure  gra¬ 
dients,  produced  by  gradients  of  liquid  saturation  and  surface 
tension  [5,7],  and  so,  convection  is  neglected  in  porous  medi¬ 
ums.  The  capillary  pressure  is  semi-empirically  calculated  by 
means  of  the  Leverett  function  J(S)  for  hydrophobic  unsatu¬ 
rated  porous  mediums  [2,7]: 


<7  cos  0 

l KfW 


J(S) 


The  Leverett  function  J(S )  for  hydrophobic  porous  mediums  is 
expressed  as 


7(5 )  =  1.4175  -  2.12052  +  1.26353 


where  the  reduced  saturation  factor  S  is  the  following: 


Below  the  saturation  value  S[m  (immobile  saturation),  the  water  is 
assumed  to  be  a  discontinuous  flow  pattern,  and  so,  the  capillary 
flow  is  neglected.  At  the  gas  channel,  liquid  water  is  disposed 
either  by  high  convective  air  flow  or  by  evaporation.  Then,  the 
reduced  saturation  factor  is  assumed  to  be  zero  at  the  gas  channel 
and  gas  diffusion  layer  interface  [7] . 

To  determine  the  phenomenological  behavior  of  the  liquid  in 
the  porous  medium,  Darcy’s  law  is  used: 


Q 

v Qu2Q(L)KpS 

MH20(L) 


VS  =  Aglh2o 


In  order  to  calculate  the  change  of  phase  of  water  to  establish 
the  saturation  factor  in  porous  medium,  the  condensation  kinetic 
theory  is  used  to  determine  the  volumetric  condensation  rate 
^glh2o  expressed  as 


^glh2o  =  KglMh2o 


Alg  ph2o(G)  ~  Pu2o 
VL  RT 


(10) 


The  factor  ^gl  consists  of  the  diffusion  terms  for  the  con¬ 
densation  process  [7].  The  evaporation  was  assumed  to  be  the 


The  density  of  gases  is  calculated  using  the  ideal  gas  equation. 
Dynamic  viscosity  and  thermal  conductivity  are  determined  us¬ 
ing  the  kinetic  theory  of  gases. 

The  thermal  conductivities  are  expressed  as  an  average  be¬ 
tween  conductivities  of  phases  involved  in  each  infinitesimal 
element  of  the  system: 

k  =  kSf(l  -  0  +  kGf(€)(  1  -  S)  +  kLf(€)s  (11) 

The  volumetric  thermal  capacity  is  also  regarded  as  an  average 
between  phases.  Consequently,  thermal  equilibrium  between 
phases  is  assumed: 


Cp  =  cps(l  -  6)  +  cpGe(l  -  s)  +  cpLes 


(12) 


In  order  to  consider  the  influence  of  the  porous  medium  on  the 
conductive  heat  transport,  the  following  correction  factors  (Eqs. 
13  and  14)  were  used  for  the  catalyst  layers  and  the  gas  diffusion 
layers,  respectively  [8]: 

m  =  J5  (i3) 

<»> 

The  factor  is  considered  as  0.1 1  and  a  is  0.527  for  in-plane 
diffusion  and  0.785  for  cross-plane  diffusion  [7]. 

The  mixture  diffusivity  was  approximated  by  the  mixture  rule 
[9]: 


D/.m  — 


J2i=l  xiMi  J2i=l,i^j(xi/ Dij) 


(15) 


The  binary  diffusivities  were  calculated  using  the  following  cor¬ 
relation  for  gases  at  low  density  [9] : 


_ _ 

(^c^c;)1/3(^TCi-)5/12((l/M/)  +  (1/M;))1/2 


(16) 


in  this  equation,  a  and  b  are  constant.  If  the  mixture  diffusivity 
is  calculated  for  a  porous  medium,  Eqs.  (13)  and  (14)  are  used 
to  account  for  the  effect  of  porosity  and  tortuosity  of  the  cata¬ 
lyst  layer  and  the  gas  diffusion  layer,  respectively.  The  effect  of 
the  liquid  water  saturation  is  accounted  for  by  using  a  normal¬ 
ized  function  f(s)  =  (l  —  s)15 .  Therefore,  the  final  form  of  the 
diffusivity  is 


^z‘,m-corrected  —  (s) 


(17) 
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For  the  polymer  electrolyte,  the  electro-osmotic  drag  coeffi¬ 
cient  Aforag  was  calculated  at  30  and  80  °C  using  the  empirical 
data  for  Nafion  from  Springer  et  al.  [10]  and  Neubrand  [11], 
respectively.  Diffusion  coefficient  of  water  D\  was  calculated  at 
30  and  80  °C  from  the  empirical  data  for  Nafion  of  Zawodzinski 
et  al.  [12]  and  Neubrand  [11],  respectively. 


2.3.  Boundary  conditions 


The  thermal  boundary  conditions  between  flow  and  channel’s 
walls  may  be  considered  as  either  adiabatic  or  isothermal. 

The  water  content  of  the  polymer  electrolyte  was  assumed  to 
be  highly  affected  by  external  conditions.  Therefore,  a  simple 
theory  was  considered  suitable  for  determining  the  interfacial 
mass  transfer  between  the  cathode  and  anode  porous  mediums 
and  polymer  membrane.  The  water  transport  between  the  cata¬ 
lyst  layer  and  the  polymer  membrane  is  driven  by  the  difference 
between  the  bulk  concentrations  in  each  phase  and  the  equilib¬ 
rium  concentrations  of  water  in  the  interface.  In  this  work  the 
interfacial  mass  transfer  relation  of  the  molar  fraction  of  each 
phase  (Yi  and  Xi)  and  the  mixture  diffusivity  of  the  species  i  by 
each  phase  (A;mi  and  Dfm 2)  is  expressed  by 


Q\Dj  m\ 


Xj  -  V* 

M\  Avi 


Q2^i,  m2 


Y*  -  Yt 


M2AX2 


(18) 


The  equilibrium  data  from  Springer  et  al.  [10]  at  30  °C  and 
from  Hinatsu  et  al.  [13]  at  80  °C  were  used  to  complete  this 
interfacial  mass  transfer  model  by  relating  the  molar  fraction  of 
water  of  each  phase  at  the  interface  (Y*  and  X*). 


2.4.  Electrochemical  model 


The  electrochemical  kinetic  was  determined  by  using  the  ag¬ 
glomerate  model  [14-18].  The  agglomerate  model  states  that 
the  catalyst  layer  is  a  porous  medium  composed  by  platinum 
and  carbon  particles  mixed  together  with  polymer  electrolyte. 
This  mixture  is  defined  as  spherical  agglomerates,  which  have 
void  spaces  between  them.  In  this  work,  the  polymer  electrolyte 
is  assumed  as  Nafion,  the  water  is  assumed  to  be  produced  in 
gas  phase  and  condensation  could  take  place  in  this  medium. 
The  anode  reaction  is  assumed  to  be  dominated  by  the  cathode 
reaction  kinetic,  and  so,  the  reaction  rate  is  finally  expressed  as 

rH+  =  gg  2  - Htfflagg  COth  e  flagg  -  1)  (19) 

^agg 

with  the  Thiele  Modulus  (0Ragg)  defined  as 


^pt^superf^'oref  „  f  OIcFt]\ 

roaggCo2ref  agg6XPV  2 Rt) 


(20) 


The  Henry  constant  for  the  dissolution  of  oxygen  in  water  is  used 
to  calculate  the  concentration  at  the  surface  of  the  agglomerate 
[19]: 


P  (666 
Co2s  -  xq2  1Q1325  exp  [  — 


(21) 


The  reference  exchange  current  density  for  the  oxygen  reduction 
in  Nafion  (/oref)  is  calculated  by  the  correlation  presented  by 
Wang  et  al.  [14]  made  from  experimental  data  [20]: 

4001 

log (f oref)  =  7.507  -  —  (22) 

The  diffusivity  of  oxygen  in  Nafion  is  calculated  from  the  cor¬ 
relation  presented  by  Marr  and  Li  [16],  which  is  made  from 
experimental  data  [20].  Then,  the  oxygen  diffusivity  in  the  ag¬ 
glomerate  Dagg  is  expressed  as 

Dagg  =  t— 1.07  X  10~5 


+  9.02  x  10  6  exp 


T-  273\\ 
106.7  J ) 


(23) 


To  account  for  the  voltage  losses  which  are  induced  by  the  Nafion 
resistance  to  the  transport  of  protons  to  the  active  sites,  a  one  di¬ 
mensional  expression  of  Ohm’s  law  is  defined  to  estimate  the  in¬ 
crease  of  the  cathode  overpotential  in  the  cathode  catalyst  layer: 

^H+((l  —  6)6agg)1‘5—  =  Nn+F  (24) 

ox 

where  kH+[^_1m“  l]  is  the  proton  conductivity  in  Nafion  ex¬ 
pressed  as  [10] 

kH+  =  (0.5139A  -  0.3260) exp  M268  (+  _  I 


To  calculate  the  fuel  cell  potential,  the  cathode  overpotential 
and  the  polymer  membrane  voltage  loss  are  subtracted  from  the 
equilibrium  potential.  The  voltage  loss  produced  by  the  bipolar 
plates  and  the  contact  resistance  between  the  gas  channel  ribs 
and  gas  diffusion  layers  was  neglected: 

E  =  1.23  —  9  x  10+T-298) 

23RT  ,  i 

+  ^  Wk'-O.)  -  >>  -  T- 


2.5.  Source  terms 


The  pressure  drop  in  gas  channels  was  estimated  using  the 
Darcy  friction  factor,  and  so,  the  implementation  of  the  SIM¬ 
PLE  algorithm  [21]  was  avoided.  Therefore,  only  diffusion  is 
assumed  to  be  the  primary  mechanism  in  gas  diffusion  layers 
and  catalyst  layers  while  convection  is  neglected.  According  to 
this  assumption,  the  source  term  of  momentum  Sy  in  the  porous 
mediums  is  not  considered  in  Eq.  (1),  and  so,  the  Navier-Stokes 
equations  are  only  used  in  the  gas  channels  without  source  terms. 
The  source  term  of  mass  Sm  in  Eq.  (4)  was  also  neglected  fol¬ 
lowing  the  demonstration  made  by  Wang  and  Wang  [22] . 

Joule  heating  in  the  polymer  membrane,  heat  produced 
in  cathode  reaction  and  heat  produced  in  condensation- 
evaporation  entropy  change  are  taken  into  account  as  source 
terms  of  energy  St  in  Eq.  (2).  The  Joule  heating  Sti  in  the  mem¬ 
brane  is  calculated  by  using  the  proton  conductivity  in  Nafion 
and  the  proton  flow  as  follows: 

•SY,  =  iNf  F)1  (27) 

kH+ 
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The  heat  produced  in  the  cathode  reaction  is  determined  using 
the  reaction  heat  of  oxygen  reduction  AHvo2  and  the  reaction 
rate: 

St2  =  AHro2ro2  (28) 

And  the  condensation-evaporation  entropy  change  in  the  cath¬ 
ode  porous  mediums  is  estimated  using  the  latent  heat  of  vapor¬ 
ization  of  water  at  local  temperature  //lg  and  the  volumetric 
condensation  rate: 

St3  =  #lg^glh2o  (29) 

The  electrochemical  reaction  in  the  catalyst  layers  is  consid¬ 
ered  as  the  source  term  of  species  i  5)  in  Eq.  (3): 

Si  =  Min  (30) 

The  reaction  rate  of  every  species  is  obtained  from  the  one  ex¬ 
pressed  in  Eq.  (19)  using  the  stoichiometrical  relations  derived 
from  the  well-known  mechanism  of  reaction  in  PEM  fuel  cells. 

Finally,  the  condensation  and  evaporation  process  in  the  cath¬ 
ode  porous  medium  is  considered  as  additional  source  term  for 
the  water  gas  balance: 

Sh2o  =  ^glh2o  (31) 

It  is  important  to  mention  that  every  source  term  is  activated 
only  in  its  defined  domain. 

3.  Results  and  discussion 

In  PEM  fuel  cells,  there  are  many  parameters  that  may  affect 
their  performance  either  slightly  or  heavily.  These  parameters 
may  be  separated  into  three  sections:  the  geometrical,  the  elec¬ 
trochemical  and  the  operational.  The  geometrical  and  electro¬ 
chemical  parameters  arbitrarily  used  in  this  work  are  listed  in 
Table  1. 

The  principal  aim  is  the  study  of  the  behavior  of  the  polymer 
membrane  water  content,  the  saturation  in  the  cathode  porous 
medium  and  the  cathode  reaction  rate,  in  order  to  establish  the 
influence  of  these  factors  on  the  performance  of  fuel  cells  un¬ 
der  different  operational  conditions.  The  results  are  mostly  pre¬ 
sented  as  transversal  averages  along  the  gas  channels,  given  that 
the  three-dimensional  displays  are  not  practical  to  make  com¬ 
parisons. 


Table  1 


Geometrical  and  electrochemical  parameters 


Description 

Value 

Channel  height  (m) 

0.0004 

Channel  width  (m) 

0.001 

Channel  length  (m) 

0.10 

Channel  rib  (m) 

0.0004 

Gas  diffusion  layer  thickness  (m) 

0.0001 

Gas  diffusion  layer  porosity 

0.5 

Catalyst  layer  thickness  (m) 

0.000001 

Catalyst  layer  porosity 

0.4 

Polymer  membrane  thickness  (m) 

0.00023 

Pt  loading  (mgptcm-2) 

1 

Pt  surface  area  (cmptgPt-1) 

250000 

Agglomerate  porosity 

0.2 

Agglomerate  radii  (cm) 

0.00001 

Reaction  heat  (J  mol-1) 

460000 

3. 1.  Water  content  field  at  different  current  densities 

The  magnitude  of  the  current  density  may  affect  the  water 
content  field  of  the  polymer  membrane.  Therefore,  the  local 
current  density  and  the  local  water  content  along  the  channel 
were  computed  at  different  total  average  current  densities  and 
fully  humidified  inlet  flows  (Fig.  2(A)  and  (B),  respectively).  In 
Fig.  2(A),  it  can  be  seen  that  the  local  current  density  dimin¬ 
ishes  along  the  channel,  because  of  the  depletion  of  oxygen  by 
the  cathode  reaction.  The  local  water  content  increases  along 
the  channel,  given  the  increase  of  water  activity,  as  a  result  of 
the  production  of  water  by  the  electrochemical  reaction  in  the 
cathode  (Fig.  2(B)).  The  reaction  at  the  anode  also  contributes  to 
the  increase  of  the  water  activity  at  this  side,  due  to  the  fully  hu¬ 
midified  condition  of  the  anode  inlet  flow  and  the  consumption 
of  hydrogen  along  the  channel.  Consequently,  the  water  molar 
fraction  increases  along  the  channel  at  the  anode  side.  The  in¬ 
crease  of  water  content  becomes  slower  towards  the  end  of  the 
channel,  as  a  product  of  the  consumption  of  oxygen  along  the 
channel  and  so  the  reduction  of  the  cathode  reaction  rate.  Due 
to  the  fully  humidified  condition  of  the  inlet  flows,  the  electro- 
osmotic  drag  is  insufficient  to  dry  the  anode  side  of  the  polymer 
membrane,  even  for  high  current  densities.  Under  these  condi¬ 
tions  the  cell  voltage  is  only  affected  by  the  proton  flux  and  not 
by  the  water  content  of  the  polymer  membrane. 


Fig.  2.  (A)  Current  density  and  (B)  membrane  water  content  along  the  channel  for  different  average  current  densities. 


208 


L.  Matamoros,  D.  Briiggemann  /  Journal  of  Power  Sources  161  (2006)  203-213 


(A) 


Channel  length  [m] 


Fig.  3.  (A)  Current  density  and  (B)  membrane  water  content  along  the  channel  for  different  relative  humidities  of  anode  inlet  flow. 


3.2.  Effect  of  the  cathode  and  anode  inlet  flow  relative 
humidity 

As  mentioned  above,  once  the  inlet  flows  are  fully  humidi¬ 
fied,  the  water  content  of  the  polymer  membrane  is  favored,  and 
consequently,  the  performance  of  the  fuel  cell  is  not  affected  by 
the  hydrating  conditions  of  the  electrolyte.  However,  the  rela¬ 
tive  humidity  of  the  cathode  inlet  flow  may  have  a  considerably 
influence  on  the  performance  of  the  fuel  cell.  In  case  the  inlet 
flows  are  not  well  humidified,  voltage  losses  from  the  electrolyte 
may  increase  and  have  a  high  negative  influence  on  the  cathode 
reaction  rate,  as  a  result  of  the  decrease  of  proton  conductivity 
of  the  polymer  electrolyte,  which  depends  linearly  on  the  water 
content.  If  the  water  content  decreases,  then  the  cathode  over¬ 
potential  increases,  due  to  the  increment  of  proton  resistance, 
and  so  the  cathode  reaction  rate.  The  voltage  loss  induced  by 
the  polymer  membrane  is  also  increased,  as  product  of  the  de¬ 
hydrating  conditions. 

The  relative  humidity  of  cathode  and  anode  inlet  flows  were 
varied  at  high  current  density  operation,  in  order  to  study  the 
performance  of  the  fuel  cell  under  dehydrating  conditions.  The 
relative  humidity  of  the  cathode  inlet  flow  is  maintained  at  30% 
and  the  anode  inlet  flow  presents  relative  humidity  of  50, 70  and 
90%.  Fig.  3(A)  shows  that  the  local  current  density  diminishes 
along  the  channel,  as  a  result  of  the  decrease  of  oxygen  con¬ 
centration.  The  effect  of  the  electrolyte  on  the  current  density  is 


slightly  seen  for  the  anode  relative  humidity  of  90%,  given  that 
the  current  density  is  higher  at  the  gas  channel  inlet.  However, 
the  oxygen  is  faster  consumed,  and  so,  the  local  current  den¬ 
sity  is  less  in  comparison  to  the  50  and  70%  relative  humidity. 
The  local  water  content  increases  along  the  channel  (Fig.  3(B)), 
consequently  the  electrolyte  limitations  at  the  cathode  catalyst 
layer  decrease  as  well.  In  this  case,  the  anode  relative  humidity 
enhances  the  hydration  of  the  polymer  membrane.  Under  these 
conditions  it  can  be  seen  that  the  effect  of  the  polymer  electrolyte 
on  the  current  density  is  small  in  comparison  to  the  effect  of  the 
oxygen  mass  transfer  limitations  and  depletion. 

To  study  the  potential  effect  of  Nafion  in  the  cathode  electrode 
on  the  reaction,  the  thickness  of  the  cathode  catalyst  layer  was 
made  thicker  (by  a  factor  of  10),  maintaining  the  superficial 
platinum  load  constant,  in  order  to  increase  the  proton  transport 
resistance  in  comparison  to  the  reaction  rate.  The  results  are 
showed  in  Fig.  4(A)  and  (B)  at  the  same  humidity  conditions 
used  for  Fig.  3.  Fig.  4(A)  shows  the  current  density  along  the 
channel  for  dehydrating  conditions.  It  can  be  seen  that  for  low 
water  content  of  polymer,  the  resistance  of  the  Nafion  to  the 
proton  transport  becomes  high  and  affects  the  performance  of 
the  fuel  cell  so  as  the  oxygen  mass  transfer  limitations  do.  In  this 
case,  the  current  density  slightly  increases  along  the  channel, 
given  that  the  water  content  of  the  electrolyte  becomes  higher 
as  result  of  the  production  of  water  by  the  cathode  reaction  (Fig. 
4(B)).  Then,  the  oxygen  mass  transfer  limitation  and  the  decrease 
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Fig.  4.  (A)  Current  density  and  (B)  membrane  water  content  along  the  channel  for  different  relative  humidities  of  anode  inlet  flow. 
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Fig.  5.  Maps  of  (A)  vapor  mass  fraction,  (B)  polymer  membrane  water  content  field  and  (C)  current  density. 


of  the  oxygen  concentration  reduce  the  reaction  rate,  despite  the 
hydration  of  the  polymer  membrane  and  the  electrolyte  in  the 
cathode  catalyst  layer.  Fig.  5  shows  the  cathode  water  vapor 
mass  fraction,  the  water  content  field  for  relative  humidity  of 
cathode  and  anode  inlet  flows  at  30  and  50%  (respectively)  and 
the  current  density.  As  it  can  be  seen,  the  maximal  current  density 
area  has  moved  towards  the  end  of  the  channel,  as  a  result  of 
the  better  hydration  of  the  electrolyte  by  the  electrochemical 
reaction.  However,  this  area  is  closed,  because  of  the  decrease 
of  the  oxygen  concentration  along  the  channel. 

3.3.  Water  saturation  of  the  cathode  porous  medium 

The  cathode  saturation  represents  one  of  the  most  impor¬ 
tant  problems  to  consider,  given  that  the  liquid  water  may  block 
the  path  of  reactants  to  the  active  sites.  Therefore,  the  inlet 
relative  humidity  improves  the  electrolyte  performance  of  the 
fuel  cell,  but  may  decrease  the  flux  for  the  diffusion  to  ac¬ 
tive  sites.  The  saturation  should  be  quantitatively  estimated  to 
avoid  the  flooding  in  the  cathode  side  by  an  incorrect  water 
management. 

The  saturation  in  the  cathode  porous  medium  may  be  affected 
by  several  factors:  current  density,  temperature  of  operation,  hy- 
drophobicity  and  permeability  of  porous  medium  and  relative 
humidity  of  cathode  and  anode  inlet  flows.  The  condensation 
and  evaporation  process  depends  highly  on  the  interfacial  area 
between  the  gas  and  the  liquid.  To  predict  an  accurate  value  for 
the  volumetric  interfacial  area  is  a  difficult  task,  especially  for  the 


complexity  of  the  three-phase  system.  As  mentioned  above,  the 
volumetric  interfacial  area  is  arbitrarily  considered  as  a  constant 
parameter  during  the  condensation  and  evaporation  process. 

The  saturation  of  the  porous  medium  becomes  a  problem 
when  the  system  is  unable  to  dispose  the  liquid  water  product 
from  the  excess  of  water  vapor  at  the  same  speed  it  is  produced 
by  the  cathode  reaction.  Therefore,  the  liquid  water  accumu¬ 
lates  and  the  available  area  for  the  diffusion  is  reduced.  Then, 
high  voltage  losses  occur,  because  of  the  high  mass  transfer  con¬ 
trol  from  the  high  mechanical  resistance  for  diffusion.  However, 
once  the  water  accumulates,  some  continuous  liquid  water  paths 
will  be  formed  and  there  will  be  a  capillary  flow  straight  to  the 
gas  channel,  so  that  the  liquid  water  can  be  disposed.  This  phe¬ 
nomenon  may  keep  the  system  from  flooding,  although  the  mass 
transfer  control  may  remain. 

In  this  work,  the  influence  of  the  volumetric  interfacial  area 
on  the  performance  of  the  fuel  cell  was  computed  by  assuming 
different  values  of  this  parameter.  Fig.  6(A)  and  (B)  show  the 
local  saturation  and  water  content  along  the  channel  for  volumet¬ 
ric  interfacial  areas  of  1000,  10,000  and  50,000  m2  m~3  at  high 
current  densities  and  fully  humidified  inlet  flows,  respectively.  It 
can  be  observed  from  Fig.  6(A)  that  the  amount  of  liquid  water  is 
slightly  governed  by  the  effect  of  the  interfacial  area,  for  the  val¬ 
ues  chosen  in  this  work.  According  to  these  results,  the  system 
is  able  to  get  rid  of  the  water  and  the  maximal  saturation  is  ap¬ 
proximately  below  the  value  of  0.2  at  the  cathode  catalyst  layer. 
Consequently,  the  capillary  flow  controls  the  saturation  process 
against  the  condensation  and  evaporation  process.  In  this  case 


Fig.  6.  (A)  Membrane  water  content  and  (B)  saturation  along  the  channel  at  different  volumetric  gas-liquid  interfacial  areas. 
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Fig.  7.  Cell  voltage  and  power  vs.  average  current  density  at  different  volumetric 
gas-liquid  interfacial  areas. 


the  dominating  effect  is  neither  the  volumetric  interfacial  area 
nor  the  current  density  of  the  fuel  cell,  but  the  permeability  and 
hydrophobicity  of  the  porous  medium.  It  can  also  be  seen  from 
Fig.  6(A)  how  the  average  saturation  changes  along  the  channel. 
This  factor  increases  as  a  result  of  the  increase  of  water  activity 
along  the  channel.  The  saturation  tends  to  reach  a  maximum, 
given  the  decrease  of  reaction  rate  towards  the  channel.  The  wa¬ 
ter  content  could  be  affected  by  the  volumetric  interfacial  area 
of  the  gas  and  liquid  interface,  as  it  can  be  seen  in  Fig.  6(B).  The 
local  water  content  decreases  with  increasing  volumetric  inter¬ 
facial  area.  Therefore,  the  higher  condensation  of  water  vapor 
reduces  the  mass  of  water  that  may  be  transferred  to  the  poly¬ 
mer  electrolyte.  Consequently,  the  voltage  loss  in  the  polymer 
membrane  may  be  affected  by  the  higher  volumetric  interfacial 
area. 

Fig.  7  shows  the  cell  voltage  and  power  vs.  average  current 
density  for  volumetric  interfacial  areas  of  1000,  10,000  and 
50,000  m2m-3.  As  it  can  be  seen  for  high  current  densities, 
fuel  cell  voltage  decreases  while  volumetric  interfacial  area  in¬ 
creases.  Therefore,  under  these  conditions  the  current  density  is 
not  seriously  affected  by  saturation  of  liquid  water  in  the  cathode 
electrode,  given  that  the  difference  of  the  saturation  magnitudes 
is  small  for  all  cases.  The  significance  of  this  behavior  may  be  the 
low  mass  transfer  control  that  liquid  water  induces  in  the  cathode 
porous  medium  in  comparison  to  the  resistance  that  the  porous 
medium  exerts  on  the  diffusive  flux.  The  performance  of  the  fuel 
cell  is  mainly  decreased  by  the  influence  of  the  volumetric  inter¬ 
facial  area  on  the  water  uptake  of  the  polymer  membrane,  and  so 
on  voltage  of  the  fuel  cell.  Nevertheless,  the  water  uptake  from 
liquid  water  is  neglected  and  this  factor  might  be  important  to 
obtain  more  accurate  results  for  the  water  uptake  of  the  polymer 
membrane. 

To  study  the  influence  of  the  cathode  inlet  relative  humidity 
on  the  liquid  saturation  of  the  cathode  porous  medium,  different 
curves  were  computed  at  30,  60  and  90%  cathode  relative  hu¬ 
midity  maintaining  the  anode  relative  humidity  at  100%  and  the 
volumetric  interfacial  area  at  1000  m2  m-3  (Fig.  8).  In  this  case 
the  relative  humidity  has  a  strong  influence  on  the  condensation 
and  evaporation  process,  due  to  the  unsaturation  of  the  inlet  flow. 


Fig.  8.  Saturation  along  the  channel  at  different  cathode  inlet  flow  relative  hu¬ 
midity. 


Then,  there  will  be  a  tendency  to  evaporation  at  the  gas  channel 
and  cathode  gas  diffusion  layer,  and  so  the  dispose  of  water  is 
made  by  change  of  phase  and  not  by  the  drag  of  liquid  droplet 
in  the  gas  channel.  For  higher  inlet  humidity  the  capillary  flow 
dominates  against  the  evaporation,  and  so  the  disposal  is  made 
by  the  mechanical  energy  of  the  flow  in  the  gas  channel. 

The  relation  of  local  temperature  and  liquid  saturation  of 
the  cathode  porous  medium  is  also  to  be  accounted  for.  Fig. 
9  shows  the  distributions  of  temperature  and  liquid  saturation 
in  the  through-plane  direction  along  the  channel  right  in  the 
middle  of  the  fuel  cell  for  fully  humidified  inlet  flows  and  high 
current  density  operation  (1.04  A  cm-2).  For  fully  humidified 
inlet  flows,  reaction  rate  is  higher  at  the  entrance  and  diminishes 
along  the  channel  by  depletion  of  reactants.  Likewise  temper¬ 
ature  is  higher  at  the  entrance  and  becomes  lower  along  the 
channel  (Fig.  9(A)).  In  through-plane  direction,  temperature  di¬ 
minishes  towards  the  gas  channel  due  to  convection  in  the  in¬ 
terface  between  gas  channel  and  gas  diffusion  layer.  According 
to  the  distribution  of  temperature,  the  convection  in  gas  chan¬ 
nels  maintains  the  cell  almost  at  uniform  temperature.  Under 
the  conditions  listed  in  Table  1,  the  heat  of  reaction  and  Joule 
heating  are  not  high  enough  to  increase  considerably  the  tem¬ 
perature  against  convection  in  gas  channels.  In  the  case  of  the 
local  saturation,  the  distribution  (Fig.  9(B))  depends  not  only  on 
temperature  distribution,  but  also  on  water  vapor  mass  fraction 
distribution  and  capillary  flow.  Fig.  9(B)  shows  that  the  local  liq¬ 
uid  saturation  slightly  increases  along  the  channel  and  slightly 
decreases  in  through-plane  direction  towards  the  gas  channel. 
Thus,  as  mentioned  above  capillary  flow  plays  the  most  impor¬ 
tant  role  on  the  liquid  distribution.  Even  though  water  activity 
may  change  almost  100%  along  the  channel  (as  it  can  be  seen  in 
Fig.  5(A)),  local  saturation  is  almost  homogeneous  in  cathode 
porous  medium  due  to  the  capillary  flow.  The  lower  liquid  satu¬ 
ration  after  the  entrance  of  fuel  cell  shows  that  temperature  and 
water  activity  distributions  have  an  influence  on  the  distribution 
of  saturation;  however  these  effects  are  small  in  comparison  to 
the  capillary  flow. 
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Fig.  9.  Through-plane  and  along-channel  distributions  of  (A)  temperature  and 
(B)  liquid  saturation. 


3.4.  Effect  of  the  operational  temperature 

To  observe  the  influence  of  the  operational  temperature  on 
the  performance  of  the  fuel  cell,  the  local  water  content  and  the 
current  density  along  the  channel  were  computed  for  70,  80  and 
90  °C  (Fig.  10(A)  and  (B),  respectively).  From  Fig.  10(A),  it  can 
be  seen  that  the  local  current  density  decreases  while  the  tem¬ 
perature  increases.  This  effect  is  especially  seen  for  the  results 
of  90  °C.  Therefore,  the  performance  of  the  fuel  cell  diminishes 
at  higher  temperatures  for  fully  humidified  inlet  flows,  as  a  re¬ 
sult  of  the  influence  of  the  temperature  on  the  hydration  of  the 
polymer  electrolyte,  as  it  can  be  seen  in  Fig.  10(B).  For  the  tem¬ 
perature  of  90  °C,  the  water  content  of  the  polymer  membrane 
dramatically  diminishes  for  fully  humidified  inlet  flows,  and  so 
the  electrolyte  overpotential  increases.  Hence,  the  temperature 
may  have  a  positive  effect  on  the  reduction  of  the  oxygen,  but 
the  higher  proton  resistance  produced  by  the  electrolyte  dehy¬ 
dration  at  high  temperatures  has  a  dominant  influence  on  the 
cathode  reduction  reaction. 


Fig.  11.  Current  density  along  the  channel  at  different  cathode  inlet  volumetric 
flows. 

3.5.  Effect  of  cathode  inlet  volumetric  flow 

At  high  current  densities,  the  demand  of  oxygen  is  high,  and 
so  the  inlet  flow  should  be  able  to  supply  the  needed  amount 
of  oxygen.  In  the  case  the  oxygen  is  gradually  depleted,  the 
performance  of  the  fuel  cell  diminishes,  due  to  the  high  cathode 
overpotential  induced  by  the  lack  of  oxygen.  Fig.  1 1  shows  the 
local  current  densities  along  the  channel  at  different  cathode 
inlet  volumetric  flows.  The  results  show  the  high  influence  of 
the  oxygen  concentration  on  the  cathode  reaction  rate.  Even  for 
the  highest  inlet  volumetric  flow,  the  drop  of  the  current  density 
along  the  channel  is  noticeable  along  the  channel.  Nevertheless, 
the  channel  length  was  arbitrarily  chosen  and  the  aim  of  this 
calculation  is  to  show  how  the  depletion  of  oxygen  may  affect 
the  performance  of  the  fuel  cell  at  different  current  densities. 

3.6.  Comparisons  to  the  experimental  data 

The  validation  of  the  data  is  an  obliged  step  of  every  either 
experimental  or  numerical  work.  Considering  that  each  author 
has  made  emphasis  on  different  aspects  of  the  PEM  fuel  cells,  the 
comparison  between  numerical  models  becomes  useless.  Given 
that  the  properties  fields  are  not  experimentally  measured  in 
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Fig.  10.  (A)  Membrane  water  content  and  (B)  saturation  along  the  channel  at  different  operation  temperature. 
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Fig.  12.  Experimental  [23]  and  numerical  polarization  curves. 


the  literature,  but  numerically  calculated,  only  the  experimental 
polarizations  (cell  voltage  versus  average  current  density)  could 
be  taken  as  reliable  sources  to  make  comparisons.  The  shape  of 
a  polarization  curve  is  highly  influenced  at  low  and  moderates 
current  densities  by  the  electrochemical  reaction  rate,  but  at  high 
current  densities,  the  voltage  loss  in  the  polymer  membrane  also 
plays  an  important  role.  These  factors  are  highly  sensible  to  the 
electrochemical  parameters  (listed  in  Table  1)  and  the  proton 
conductivity  and  thickness  of  the  polymer  membrane,  and  so, 
it  is  necessary  to  point  out  that  the  success  of  the  validation 
depends  on  the  use  of  the  correct  electrochemical  input  data  and 
the  correct  polymer  membrane  characteristics.  Additionally,  the 
geometrical  parameters  should  be  as  close  as  possible  to  the 
parameters  of  the  experimental  setup  used  as  reference.  In  this 
work,  the  numerical  input  data  were  adjusted  as  close  as  possible 
to  the  conditions  of  the  experimental  data  from  [23],  in  order  to 
validate  the  numerical  solution.  The  comparison  between  the 
experimental  polarization  curve  at  80  °C  and  200  kPa  and  the 
numerical  one  (Fig.  12)  shows  that  the  numerical  solution  is 
suitable  for  studying  the  water  and  heat  behavior  in  PEM  fuel 
cells. 

4.  Conclusions 

A  three-dimensional  model  for  a  PEM  fuel  cell  was  pro¬ 
grammed  in  Fortran90  to  achieve  a  numerical  solution  of  the 
water  and  heat  behavior,  in  order  to  analyze  the  factors  that  may 
be  either  positive  or  negative  for  the  performance  of  this  tech¬ 
nology.  Under  the  conditions  and  geometry  chosen  in  this  work, 
the  system  showed  to  be  highly  dependant  on  the  hydration  of 
the  polymer  membrane,  which  was  affected  by  the  relative  hu¬ 
midity  of  the  inlet  flows.  The  saturation  of  the  cathode  porous 
medium  was  dominated  by  the  capillary  flow,  and  so,  the  per¬ 
meability  was  high  enough  to  dispose  the  water  in  the  system. 
The  volumetric  interfacial  area  of  the  water  affected  the  perfor¬ 
mance  of  the  fuel  cell  by  reducing  the  amount  of  water  to  be 
transferred  to  the  polymer  membrane.  The  operational  tempera¬ 
ture  may  have  a  negative  effect  on  the  hydrating  properties  of  the 


polymer  electrolyte,  and  so  on  the  performance  of  the  fuel  cell, 
especially  for  values  above  90  °C.  The  saturation  of  the  cathode 
porous  medium  was  not  dramatically  affected  by  the  operation 
conditions  and  its  magnitude  did  not  overcome  the  value  of  0.2. 
Then,  the  capillary  flow  was  high  enough  to  dispose  the  produced 
water  and  to  prevent  the  system  from  flooding.  The  chosen  ex¬ 
perimental  polarization  curve  was  in  good  agreement  with  the 
calculated  one.  This  numerical  solution  constitutes  a  useful  tool 
for  studying  the  behavior  of  PEM  fuel  cell  under  several  condi¬ 
tions.  The  characteristics  of  the  materials  that  compose  “ME A” 
can  also  be  varied  to  extend  the  possibilities  of  this  solution.  For 
further  research,  this  numerical  solution  can  be  adapted  to  more 
complex  geometries  (serpentine  and  interdigitated  flow  fields). 
The  estimation  of  the  efficiency  for  different  configurations  is 
also  an  aspect  to  consider  for  future  aims. 
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